Cloud forest plant assemblages are among the richness worldwide. Hummingbird visited plants make up a large percentage of the plant assemblage. Previous work has shown that hummingbirds feed on plants that match corollas based on bill length. This suggests that plants with similiar corolla morphologies compete for pollination resources. To minimize competition between species, as well as decrease amount of pollen loss due to heterospecific transfer, there should be selection for minimal overlap in flowering timing. While this has been evaluated within a single genus (Stiles 1975) there has been little work in combining relatedness, environment and morphology in a single framework. My goal is to compare the effect of morphology of competiting species on the abuncance of flowers along a wide elevation gradient.
We collected flower abundance for hummingbird visited flowers along a 1300m elevation gradient in northern ecuador.
We measured three flower morphology traits for each species: corolla width, corolla length and height of plant from ground.
Flower abundance changes throughout the year and along the elevation gradient.
\[Abundance \sim Elevation + Julian Day + Error\]
\[Error \sim \text {Morphology of Co-occurring Species}\]
The goal of the analysis is to measure the effect of this error term on shaping abundance of species flowers. Additional analysis could also look at elevation/distance to a conspecific as promoting abundance.
The abundance of plants is models as a log normal distribution instead of a poisson to have control over the variance term.
\[Abundance \sim LogN(\mu,\sigma) \]
This analysis needs to broken down into several tractable steps.
\[\lambda \sim Intercept + \beta_1 * Elevation \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \]
\[\mu \sim Intercept + \beta_1 * Elevation \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[\sigma \sim Gamma(0.0001,.0001) \]
Sigma is drawn as a gamma, because the conjugate of lognormal is a gamma.Hopefully this will help convergence.
\[Abundance \sim LogN(\mu,\sigma) \]
\[\mu \sim Intercept + \beta_1 * Elevation + \beta_2 * \text{Julian Day}. \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[\sigma \sim Gamma(0.0001,.0001) \] \[ \text{Julian Day} \sim U(0,365) \] \[\beta_2 \sim N(0,.0001) \]
\[Abundance \sim LogN(\mu,\sigma) \]
\[\mu \sim Intercept + \beta_1 * Elevation + \beta_2 * \text{Julian Day}. \]
\[intercept \sim N(0,.0001) \] \[\beta_1 \sim N(0,.0001) \] \[ \text{Julian Day} \sim U(0,365) \] \[\beta_2 \sim N(0,.0001) \]
With the errors: \[\sigma \sim Gamma(lambda,.001)\] \[ \lambda \sim Intercept + \beta_3 * \text{Morphology of Co-occurring Species}\] \[ \beta_3 \sim N(0,0.0001)\]
d<-read.csv("C:/Users/Ben/Dropbox/Thesis/Maquipucuna_SantaLucia/Results/FlowerTransects/FlowerTransectClean.csv",row.names=1)
#remove missing or infinite data
d<-d[is.finite(d$Total_Flowers),]
d<-d[is.finite(d$ele),]
Sampling structure throughout the year
ggplot(d,aes(y=ele,x=as.Date(d$Date_F))) + geom_point() + theme_bw() + scale_x_date() + labs(x="Date")
ggplot(d,aes(x=ele,y=Total_Flowers)) + geom_point() + geom_smooth(method="lm") + labs(x="Elevation",y="Flower Abundance")
#############################
# Step #1: Load all the libraries
#############################
library(R2jags)
library(abind)
#############################
# Step #2: Set the working directory, and the directory where JAGS lives
source("Model1.R")
#############################
# Step #4: Make a list where you include all the data the model will need to run
#############################
Dat <- list(
y=as.integer(round(d$Total_Flowers,2)),
Elevation=d$ele
)
#############################
# Step #5: Make a function (with no inputs) where you put a list of parameters and their initial values
#############################
InitStage <- function() {list(alpha=1500,beta=0.1)}
#############################
#Make a column vector with the names of the parameters you want to track
#############################
ParsStage <- c("alpha","beta")
#############################
#Set the variables for the MCMC
#############################
ni <- 10000 # number of draws from the posterior
nt <- 1 #thinning rate
nb <- 2000 # number to discard for burn-in
nc <- 2 # number of chains
#############################
#Run the jags function to run the code
#############################
m = jags(inits=InitStage,
n.chains=nc,
model.file="flower.jags",
working.directory=getwd(),
data=Dat,
parameters.to.save=ParsStage,
n.thin=nt,
n.iter=ni,
n.burnin=nb,
DIC=T)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 16414
##
## Initializing model
m
## Inference for Bugs model at "flower.jags", fit using jags,
## 2 chains, each with 10000 iterations (first 2000 discarded)
## n.sims = 16000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75%
## alpha 1.459e+03 10.719 1.440e+03 1.449e+03 1.458e+03 1.468e+03
## beta -5.880e-01 0.004 -5.950e-01 -5.920e-01 -5.880e-01 -5.840e-01
## deviance 4.915e+06 38627.558 4.849e+06 4.882e+06 4.915e+06 4.948e+06
## 97.5% Rhat n.eff
## alpha 1.477e+03 1.068 28
## beta -5.810e-01 1.068 28
## deviance 4.980e+06 1.068 28
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 718434640.8 and DIC = 723349804.4
## DIC is an estimate of expected predictive error (lower deviance is better).
require(reshape2)
pars<-melt(m$BUGSoutput$sims.array[,,c('alpha','beta')])
colnames(pars)<-c("Draw","Chain","parameter","estimate")
require(ggplot2)
ggplot(pars,aes(x=Draw,y=estimate,col=factor(Chain))) + geom_point(size=.5) + geom_line() + ggtitle("Chains") + facet_wrap(~parameter,scales="free") + theme_bw() + labs(col="Chain")
require(ggplot2)
ggplot(pars,aes(x=estimate)) + geom_histogram() + ggtitle("Estimate of parameters") + facet_wrap(~parameter,scale="free") + theme_bw()